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Strong and weak thermalization of infinite non-integrable quantum systems 
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When a non-integrable system evolves out of equilibrium for a long time, local observables are 
expected to attain stationary expectation values, independent of the details of the initial state. 
However, intriguing experimental results with ultracold gases [l], [3 have shown no thermalization 
in non-integrable settings, triggering an intense theoretical effort to decide the question 3-8]. Here 
we show that the phenomenology of thermalization in a quantum system is much richer than its 
classical counterpart. Using a new numerical technique, we identify two distinct thermalization 
regimes, strong and weak, occurring for different initial states. Strong thermalization, intrinsically 
quantum, happens when instantaneous local expectation values converge to the thermal ones. Weak 
thermalization, well-known in classical systems, happens when local expectation values converge to 
the thermal ones only after time averaging. Remarkably, we find a third group of states showing no 
thermalization, neither strong nor weak, to the time scales one can reliably simulate. 
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In the 19th century, statistical mechanics was devel- 
oped as a microscopic model explaining the fundamental 
results of thermodynamics. Starting in the early 20th 
century, quantum statistical mechanics has been devel- 
oped to describe thermodynamics of quantum systems. 
However, current experiments |l|,|2j with ultracold atoms 
have attained a level of isolation and control of param- 
eters that forces us to address not just the equilibrium 
properties of quantum systems but also the question of 
how and why these systems relax to equilibrium starting 
from a nonthermal state. For classical systems, Boltz- 
mann's molecular chaos assumption provides a quantita- 
tive tool to describe this relaxation, and ergodic prop- 
erties of the system give the explanation. If the time 
dynamics of a system explores all states with a given en- 
ergy with uniform probability, then the long-time average 
of any observed quantity will approach the expectation 
value for that quantity in the microcanonical ensemble. 

Similarly, the emergence of a thermal bath [9( in a 
quantum system has been justified under the assump- 
tion that the long-time average of typical initial states 
produces a density matrix that approaches the thermal 
average [7|, |9j, [lCf. It has been proposed, however, that 
thermalization ma y h appen without any time average in 
the quantum case 
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indeed, it is possible that, due to 
quantum entanglement, starting from fixed, non-thermal 
initial conditions, the reduced density matrix at a given 
time t on a given region A, such that A is small com- 
pared to the system size, will converge to the thermal 
expectation value at long times t. Such a phenomenon, 
which we call "strong thermalization" , cannot occur for 
a classical system as it relies on the quantum mechanical 
fact that even if the global density matrix is a pure state, 
the reduced density matrix may be a mixed state. 

For integrable systems the existence of local conserved 
quantities prevents relaxation to the thermal state, which 
is constrained only by energy. Instead, it has been sug- 



gested that these systems will relax to a state described 
by a generalized thermal ensemble, compatible with the 
set of conserved quantities [3j, |fj|, |7j . In an infinite non- 
integrable system, there exist an infinite number of con- 
served quantities, such as the powers of the Hamiltonian. 
However, not being local, they are not expected to p re- 
vent the thermalization of local observables [7|, |9|, [ll| . 

Remarkably enough, a recent experiment |l| observed 
no signs of thermalization after long evolution in a nearly 
integrable case. The results could be in part understood 
with the hypothesis of relaxation to a generalized thermal 
ensemble |3[ in an integrable case, but from the theoret- 
ical point of view it is not clear why no thermalization 
was present away from true integrability. 

Various theoretical studies have tried to elucidate the 
question of whether or under which conditions thermal- 
ization will occur in non-integrable models [J, la, [D, |8|, Il2i - 
14| | . The study has to resort to numerical techniques, 
given the lack of analytical solutions. Moreover, the 
intrinsic computational complexity of simulating large 
quantum systems limits the affordable studies to finite 
systems or short times, so that results showing non- 
thermalization cannot be extrapolated to infinite times 
or system sizes. 

Most of the studies have tried to decide whether a given 
model thermalizes or not as a function of the Hamiltonian 
parameters [J, la, la [l2| . Recent works la, [la] . have also 
analyzed the link between appearance of quantum chaos 
and thermalization in these systems. 

Here we consider the question of relaxation using a 
fixed, non-integrable Hamiltonian, but a range of ini- 
tial states. We discover a rich phenomenology of dif- 
ferent relaxation regimes. The analyzed system is an in- 
finite translationally invariant spin chain with a nearest- 
neighbor interaction. Thanks to a recently developed 
algorithm [17|, we are able to explore its dynamics at 
relatively long times. Furthermore, the method gives us 



access to the whole reduced density matrix for a block of 
few particles. Instead of analyzing the behavior of indi- 
vidual expectation values, as previous numerical studies, 
we may then quantify the degree of thcrmalization as 
the distance between the reduced density matrix of the 
evolved state and the thermal state, pth{P) = tr L ?H ) ■ 
corresponding to the same energy. We estimate this dis- 
tance both for the instantaneous reduced state, p(t), and 
for the time averaged one p(t) = -^ J \ p{r)dT J33I ] . 

We find two clearly distinct thcrmalization regimes. 
Some states show strong thermalization. Their instanta- 
neous expectation values converge (fast) to the thermal 
ones, without the need to consider the time average. For 
other states we encounter instead weak thermalization. 
We have strong numerical evidence that such states do 
not relax even at very long times. Nevertheless, the time 
averaged observables do relax to the thermal values. Fi- 
nally, we also obtain weaker numerical evidence that for 
some third group of states no relaxation occurs, at least 
to the longest time scales we can reliably simulate. Strik- 
ingly, this situation occurs even in regions of energy for 
which the spectrum of the system is shown to be chaotic 
(see Supplementary material). Since our analysis is nu- 
merical, it may be argued that the conclusions are only 
valid for a limited range of time. However, the study 



of the properties of the method [17J, as well as the de- 
tailed analysis of errors in the current study (described 
in the Supplementary material) ensure that the reliabil- 
ity of our results extends to longer times than that of 
previous studies, enough to give strong evidence about 
the existence of these different quantum thcrmalization 
regimes. 

The first regime, or strong thermalization, is illustrated 
by the initial state |y+) (Fig. [1]), in which all spins are 
aligned along the positive y direction, corresponding to 
zero initial energy, and thus /3 = 0. The density matrix 
at any time converges fast to the thermal distribution at 
the same energy. As a consequence, the time averaged 
density matrix also converges to the thermal values. 

We identify an example of the second regime, or weak 
thermalization in the initial state \Z-\-) (Fig. [2]), in which 
all spins are initially aligned along the positive z direc- 
tion. For this state, whose energy matches that of a 
thermal state with j3 « 0.7275, no relaxation at all is 
observed for long times. Instead, the distance between 
the evolved and the thermal reduced density matrices 
seems to oscillate strongly for arbitrarily long times. All 
expectation values of the form (o~ a <8> o~b <8> cr c ) show simi- 
lar non-damped irregular oscillating behavior. Since the 
entropy increases linearly with time in the \Z+) state, 
the behavior is not due to a lack of propagating excita- 
tions (see Supplementary material). The time averaged 
density matrix shows instead an only slightly oscillating 
behavior, compatible with a very slow convergence to the 
thermal state, at a rate l/\/i. This rate is characteris- 
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FIG. 1: Strong thermalization: initial state \Y+). The main 
plot shows the distance between the three-body reduced den- 
sity matrix at instant t and the corresponding thermal state 
(/? = 0). The error bars show the difference between the result 
with the largest bond dimension (D = 120) and with a lower 
one (D = 60), which gives us an estimate of the truncation 
error. The right inset shows the distance between the time- 
averaged reduced density matrix and the thermal state. We 
superimpose a fit to a curve decaying like b/yt for long times 
(solid black line) . The left inset shows the difference between 
the thermal expectation values and the time dependent single 
body observables (a x ) (blue solid line), (a v ) (dashed green) 
and (cr z ) (dash-dotted black line). Convergence to the ther- 
mal state is observed in all three plots. 



tic of diffusive relaxation in one dimension, and has been 
seen elsewhere [12 [ . 

Finally, if the system is started in the \X+) state 
(Fig. [3|), with all spins aligned along the positive x di- 



rection, positive energy per particle, and corresponding 
to a Gibbs state with (3 sa —0.7180, the reduced den- 
sity matrix shows initially a fast relaxation. However at 
long times the distance between the evolved state and 
the thermal one is different from zero, signaling that one 
or more expectation values of local observables have not 
converged to the thermal averages. We observe that the 
density matrix has not reached a stationary value, either. 
Given the numerical character of the study, this obser- 
vation could indicate an absence of thermalization, but 
also a much slower one. It is an intriguing case, worthy 
of further, ideally experimental, study. Even in average, 
we detect no thermalization of the \X+) state, up to the 
time scales we are able to explore. 

The appearance of the different thcrmalization regimes 
does not require any fine tuning of the initial conditions 
or Hamiltonian parameters. On the contrary, by analyz- 
ing the evolution of different sets of product initial states, 
rotating the initial polarization of the spins from y to 
z, we see both strong and weak thermalization regimes 
over a range of parameters, separated by a transition. 
As shown in Fig. [H we observe that weak thermaliza- 
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FIG. 2: Weak thermalization: initial state \Z+). The dis- 
tance between the reduced density matrix for three sites and 
the thermal state of the same energy (/3 = 0.7275) oscil- 
lates strongly with time. The right inset shows the distance 
between the time averaged reduced density matrix and the 
thermal ensemble. The superimposed solid line, fit to a curve 
which behaves as b/y/i for long times, shows that the behavior 
is compatible with a slow convergence only in average. The 
left inset shows the oscillations of the individual single-body 
time dependent expectation values around the thermal ones. 
No sign of damping is observed for the longest times (t ~ 18) 
we have simulated. These results were obtained with bond 
dimension D = 240, and the error bars show the difference to 
D = 120 results. 



tion appears approximately half-way (corresponding to 
/? « 0.3502), with strong thermalization for those states 
closer to \Y+) . Similarly, we observe a transition between 
the strong and the non-thermalizing cases (see Supple- 
mentary material). 

The different regimes survive also under changes to 
the Hamiltonian parameters, showing that they are also 
not an isolated phenomenon for the particular values we 
chose. We have checked (see Supplementary material) 
that weak thermalization becomes prevalent if we de- 
crease the magnitude of the transverse magnetic field g, 
while as we increase it, strong thermalization shows up 
in a larger fraction of the initial states. The same effect 
occurs when we decrease the magnitude of the parallel 
magnetic field, h. 

Our data show that the thermalization process of a 
quantum system is a much richer phenomenon than its 
classical analogue. The appearance of the various ther- 
malization regimes cannot be linked exclusively to the 
intcgrability of the Hamiltonian, since we observe these 
different regimes for a fixed model, not close to any inte- 
grable limit. In particular, some initial states thermalize 
strongly, i.e. at the level of the instantaneous expecta- 
tion values, while others do it weakly, or only in average. 
Some other states seem to retain memory of the initial 
configuration for much longer and no type of relaxation 
could be proved. The non-relaxing configurations would 



FIG. 3: No thermalization observed: initial state \X+). The 
plot shows the time dependence of the distance between the 
evolved reduced density matrix for three sites and the ther- 
mal state (/3 = —0.7180). No thermalization is observed in 
the instantaneous density matrix, nor in the time averaged 
one (right inset). On the latter, we superimpose a fit of the 
computed values to a time dependent function, which asymp- 
totically tends to a constant (0.03). The plotted results cor- 
respond to a bond dimension D — 240, while the distance to 
the results with D = 120 is shown as error bars. In the left 
inset we plot, for the observables that determine the one site 
reduced density matrix, (<j x ) (blue solid line), (a y ) (dashed 
green) and (a z ) (dash-dotted black), the difference with re- 
spect to the thermal values. We observe that (a x ) is the 
responsible for the lack of thermalization, while all the other 
expectation values converge to the thermal averages. Study- 
ing the evolution of only a few local expectation values may 
not suffice then to detect the nonthermalization. 
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FIG. 4: Distance between the averaged 3-body reduced den- 
sity matrix and the thermal state as a function of time for var- 
ious product initial states, from \Z+) (f3 = 0.7275) to \Y+) 
(ft — 0),as indicated on the Bloch sphere on the left. 



be good candidates for an experimental study of ther- 
malization. 



Methods 

We consider an infinite translationally invariant spin 
chain with nearest-neighbour interactions of the Ising 



type, plus a local magnetic field with transverse (g) and 
parallel (h) components to the two-body interaction. 
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With a parallel (g = 0) or transverse (h = 0) magnetic 
field, the model is exactly solvable, but at a different 
angle, we have a non-integrable model, i.e. the energy 
density is the only conserved local quantity [341 ]. Wc 
simulate numerically the time evolution of various initial 
configurations under fixed Hamiltonian parameters, g = 
— 1.05 and h = 0.5, far from any integrability limits. As 
initial configurations, we choose translationally invariant 
product states. They are determined by the state of an 
individual spin, 
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In particular, the representative states discussed above 
correspond to parameters = j, (j> = (\X+)), 9 = <p = 
| (|Y+» and 6 = (|Z+». 

Using the recently developed folding method [13], we 
compute all the time dependent expectation values of 
one-, two- and three-body operators, for each initial 
state, what allows us to reconstruct the whole reduced 
density matrix for up to three sites. The thermal state 
Pth(P) with the same energy as the initial state is also 
calculated numerically with the same method (see Sup- 
plementary Material) . The distance between the evolved 
and thermal reduced density operators is then measured 
by the operator norm of their difference, d(p\,p2) = 
1 1 Pi — P2W0P1 which in this case coincides with the maxi- 
mum eigenvalue in absolute value of the difference p\ — P2 ■ 
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FIG. 5: Level spacing distribution of the unfolded spec- 
trum [TJ| for the Hamiltonian [2] in a finite (L — 14) sys- 
tem, for different Hamiltonian parameters, the superimposed 
curves show a Poissonian (dashed green) and Wigner (solid 
red) distribution, characterizing integrable and chaotic sys- 
tems, respectively. 



SUPPLEMENTARY MATERIAL 
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FIG. 6: Non integrability for a weak thermalizing state. The 
left plot shows the level spacing distribution, in the finite 
system of length L = 14, for a window of energy around 
E/N = -0.93, of width AE/N = 0.21. We observe that 
the statistics is closer to Wigner-Dyson (black) than Poisson 
(green dashed). Nevertheless, the product state with the same 
energy per site in the infinite chain shows weak thermaliza- 
tion. The right plot illustrates this, by showing the distance 
between the thermal reduced density matrix and the time 
averaged one as a function of time. The magenta diamonds 
show the data for this particular state. The strongly thermal- 
izing state \Y+) (black triangles) and the weakly thermalizing 
\Z+) (red crosses) are also shown for reference. 



THE HAMILTONIAN 



The model we consider is an infinite translationally 
invariant spin chain with an Ising type nearest neighbour 
interaction, plus a magnetic field 
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When 5 = 0, the Hamiltonian is trivially solvable, while 
for h — it reduces to the integrable Ising model with 
a transverse field, for which the exact solution can be 
found by fcrmionization. In any other case, the model is 
non-intcgrablc. 

We fix the initial values of the Hamiltonian parame- 
ters for our study to g = —1.05 and h = 0.5, which is not 
close to either of the integrable situations. In this situ- 
ation, the energy density is the only conserved quantity. 
We have additionally analyzed the spectral properties of 
the Hamiltonian, for a finite system of length L, to check 
the non-integrability in the sense of a spectrum with the 
characteristics of a random matrix ensemble. As shown 
in Fig.[5j already for L = 14 the level spacing distribution 
evidences the non-integrability of the chosen Hamiltonian 
also from the point of view of its spectral properties. For 
comparison, the level spacing distributions in both inte- 
grability limits are also shown. 

It could be argued that the non thermalization we ob- 
serve occurs for states which lie close to the edges of 
the spectrum, as \Z+) and |X+), and at these energies 
there could be an integrable effective model [20( , even for 
g = —1.05, h = 0.5. The spectral properties discussed in 
Fig.[5]would be dominated by the central part of the spec- 
trum and not reflect the properties at very low or high 
energies, while the spectrum in an energy interval in these 
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FIG. 7: Non integrability for a non-thermalizing state. The 
left plot shows that the level spacing distribution, in the finite 
system of length L = 14, for a window of energy centered 
on E/N = 0.91 and of width AE/N = 0.21, is closer to 
Wigner-Dyson than Poisson. Nevertheless, the product state 
with the same energy per site in the infinite chain shows no 
thermalization, as the right plot illustrates. The magenta 
diamonds represent the distance between the thermal reduced 
density matrix and the time averaged p for this particular 
state as a function of time. The strongly thermalizing state 
\Y+) (black triangles) and the non thermalizing \X+) (red 
crosses) are also shown for comparison. 



regions should show a very different behavior. To discard 
the integrability of the system in the interesting cases, we 
have checked the level statistics of a small energy window 
around the energy per site of a weak thermalizing state 
(Fig. [6]) and a non-thermalizing one (Fig. [7]) for the case 
of a finite system (L = 14) which can be exactly solved. 
In both cases we have found that the level spacing distri- 
bution is typical of a non-integrable system, also in these 
regions of energy. 
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FIG. 8: In the folding approach (b), operators for the same 
time step are grouped together in a double effective opera- 
tor [H. 



THE NUMERICAL METHOD 

The time evolution of an infinite ID quantum sys- 
tem is simula ted numerically within the Matrix Product 
(MPS) formalism, using the new technique 
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introduced in [17| . With this folding method, it is possi- 
ble to study the out-of-equilibrium dynamics of the sys- 
tem after longer times than with other similar methods. 

In this method, any time dependent expectation value 
(\I/(t)|0|\I/(t)) can then be expressed as a two dimen- 
sional tensor network (Fig. 8(a) ), constructed from the a 
Suzuki- Trotter expansion of the evolution operator [27 . 
|28j. Each discrete time step corresponds then to a se- 
quence of matrix product operators (MPO) [29j. In 
contrast to the standard MPS algorithm, in which the 
evolved state is approximated by an MPS after each time 
step by means of successive truncations 30, |3l( , the new 
algorithm performs the contraction of the tensor network 
in the transverse direction, i.e. along space. The left and 
tight semi-infinite lattices can then be effectively sub- 
stituted by the left and right dominant eigenvectors of 
the transfer matrix of the evolved state, (L\ and \R). 
Before contracting we apply a folding to the network 
along the time direction (Fig. |8(b)[ ). The folding oper- 
ation can be understood as performing the contraction 
(*(i)|0|*(i)) = ($|(0|*(i)) <g> |*(*)», where |$) is a 
product of unnormalized maximally entangled pairs be- 
tween each site of the chain and its conjugate (see 17| for 
a detailed discussion of the algorithm) . This is equivalent 
to grouping together tensors that correspond to the same 
time step in V? and its Hermitian conjugate, and achieves 
a more efficient representation of the entanglement in the 
transverse direction, which in turn gives access to the 
simulation of longer evolution times. 

The folding technique is appropriate for the simulation 
of time evolution, but using imaginary time evolution, it 
is also possible to efficiently compute any local expec- 
tation value in a thermal state. We obtain in this way 
the dependency of the thermal state energy with the in- 
verse temperature (3, E t h{(3) (Fig. O- From these data 
we may compute, for each one of the initial states we 
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FIG. 9: Energy per particle as a function of the inverse tem- 
perature /3 in the thermal state for the non-integrable model 
with g — —1.05, h = 0.5. The plot shows the results of the 
fourth order decomposition with 8 — 0.02 and bond dimen- 
sions D = 10 (blue crosses) and D — 20 (green crosses), and 
with step 8 — 0.01 and bond dimensions D — 20 (pink dots) 
and D — 40 (black x). We also show the results of the exact 
numerical calculation for a finite system of size L = 8 (dotted) 
and L — 12 (dashed line). 



have studied, the value of (3 corresponding to the ther- 
mal state with the same energy. This determines the 
state to which the initial configuration would be expected 
to relax, being energy density the only conserved quan- 
tity constraining the relaxation. It is then possible to 
compare the iV-particle reduced density matrix for the 
evolved state and for the thermal state at j3, to compute 
the desired distance between density operators. 

To determine the reduced density matrix for N sites, 
both for the thermal and the evolved state, we need to 
compute the expectation values of all TV-term products 
of Pauli matrices and the identity. Both in the cases of 
the thermal and the evolved state, the computation of 
the dominant left and right eigenvectors is common to 
all such expectation values. We show here the results for 
N = 3 sites. 



ACCURACY OF THE RESULTS 

Due to the numerical nature of the study, all the results 
are subject to errors. In order to assess the validity of our 
conclusions, we discuss here the character and magnitude 
of these errors using different criteria. 

The numerical method has two sources of error. The 
first one is the Trotter decomposition. The approxima- 
tion of the evolution operator by a product of exponen- 
tials introduces an error, which can be controlled by ei- 
ther reducing the parameter 8, which determines the size 
of the time step, or by using a Suzuki- Trotter expansion 
which is exact to a higher order in 8. Both decreasing 
the time step or increasing the order of the expansion 



involve longer vectors in the transverse direction, which 
will potentially worsen the truncation error. It is then 
convenient to find a trade-off between both factors. In 
our analysis we use a fourth order Suzuki- Trotter decom- 
position with time step 6 — 0.1. We have checked the 
convergence of our results with this order and time step. 

A second, generally less benign source of error is trun- 
cation, i.e. the fact of approximating a given vector by 
a MPS of fixed bond dimension, which constitutes the 
main source of numerical errors in any MPS algorithm. 
In our scheme the evolved state is not explicitly trun- 
cated, but truncation takes place in the transverse direc- 
tion, when we approximate the dominant left and right 
eigenvectors of the transfer matrix by MPS. In a stan- 
dard MPS algorithm for time evolution, truncation er- 
rors are dramatic, in the sense that, when they appear, 
the results deviate abruptly from the exact ones, and it 
becomes imposiblc to extract information from the com- 
puted quantities as soon as truncation error sets in [32J. 
However, in the folding technique, we may extract in- 
formation from longer time simulations, because errors 
come in smoothly and, even when some truncation er- 
ror occurs, the method is still able to provide a qualita- 
tive description of the physics, since we expect that our 
predictions deviate smoothly from the exact values (see 
Fig. IHil and discussion in [17|). 

To bound the numerical errors in the non-integrable 
model we may compare our results to those from other 
approaches, check the convergence of the results as we in- 
crease the bond dimension, or make use of some external 
physical criterion to assess the consistency of the com- 
puted numbers. We have used all three kinds of tests. 
First, we have cross-checked the results of our simula- 
tions with some large bond dimension simulations using 
the iTEBD algorithm [31[, in which contraction is done 
in the time direction. As shown in Fig. [lOl the fold- 
ing results with D = 120 are accurate to the longest 
times we can simulate with iTEBD, t ~ 9 — 10. Most 
remarkably, this bond dimension is enough to get a qual- 
itative description (precision 1%) of the dynamics to even 
longer times. Second, to witness the appearance of trun- 
cation errors, we have run the simulations with increas- 
ing bond dimension. The comparison of our results with 
highest bond dimension D = 240 with those for D = 120 
gives us a bound on the error, which we represent on the 
plots as an error bar. 

Finally, as a physical check of the consistency of the 
results, we test the conservation of energy along the evo- 
lution. We study the unitary evolution of a closed sys- 
tem, and the energy per particle must be constant. How- 
ever, the numerical implementation does not enforce this 
condition. On the contrary, from the Suzuki- Trotter ex- 
pansion to the truncation of the dominant eigenvectors, 
the numerical errors will in general violate this condition. 
A very large deviation between the time-evolved energy 
and the initial one would warn us about the validity of 




FIG. 10: As studied in detail in [D}], for the non-integrable 
model, and initial state JX+), magnetization ({<Tx(t))) results 
with the folded approach for D=60 (blue dots), 120 (green 
crosses), 240 (black stars) are compared to those of iTEBD 
(dash-dotted blue line for D = 256, dashed magenta for D — 
512 and black solid line for D = 1024). In the inset, we show 
the required value of D as a function of time, for different 
levels of accuracy. Notice that a moderate bond dimension 
(D < 100) suffices for a qualitative description (e ~ 1%). 



the results. As shown in Fig. QT] and [12] for the initial 
states \Z+) and |A+), with bond dimension D = 240, 
the relative error in the energy for the range of times we 
are analysing (respectively t ~ 18 and t ~ 12) is kept to 
only a few percent, consistent with the estimated trun- 
cation error. For the initial state |Y"+), with zero initial 
energy, we plot instead the expectation value of energy 
as a function of time (see Fig. [T3|) . 

One may think that this deviation of the energy could 
also introduce an error in the distance we are computing, 
as the thermal states corresponding to the computed en- 
ergy density and to the initial one will be different. To 
bound this error, we have computed the distance between 
such pair of thermal states, corresponding to the initial 
state and to the largest value of the energy found in the 
evolution, to the range of times we are showing. We find 
that this distance is significantly smaller than the one 
we observe during the dynamical evolution. In partic- 
ular, for the \Z+) initial state, the deviation in energy 
at times f w 10 reaches a 1%, which corresponds to a 
distance d(p(j3o)i p{P')) ~ 9 x 10~ 3 , while the observed 
distance between the thermal state and the evolved one 
oscillates around 0.15. For the longest times we show 
t w 18, the largest distance grows to a maximum value 
of 0.06. For the \X+) initial state, the maximum de- 
viation in energy, at t « 12, corresponds to a distance 
d(p((3o),p((3')) m 11 x 10 -3 , while the one we find at this 
same long time is 0.04 
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FIG. 11: Relative error in the energy per particle as a func- 
tion of time, for the initial state \Z+) (corresponding to 
13 = 0.7275), for bond dimension D = 40 (black crosses), 
D = 60 (blue x), D = 120 (green dots) and D = 240 (red 
stars). 
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FIG. 13: Energy per particle as a function of time, for the 
initial state \Y+) (corresponding to ft — 0), for bond dimen- 
sion D = 40 (black crosses), D = 60 (blue x) and D = 120 
(green dots). 
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FIG. 14: Entropy of the half-chain as a function of time for 
initial state \Z+) (left) and \X+) (right). The entropy is 
computed from the iTEBD simulation with different bond 
dimensions, indicated on the figures. 



Different initial states 

Our initial configurations, translationally invariant 
product states, arc specified by the state of a single spin, 



FIG. 12: Relative error in the energy per particle as a func- 
tion of time, for the initial state \X+) (corresponding to 
/3 = —0.7180), for bond dimension D = 40 (black crosses), 
D = 60 (blue x), D = 120 (green dots) and D = 240 (red 
stars) . 



DETAILED RESULTS 



Here we compile our results using various initial states 
and Hamiltonian parameters, to show the survival of the 
different thcrmalization regimes over a range of parame- 
ters. 



|tf)=cos-|0)+e^sin-|l), 

which can be represented on the Bloch sphere by the 
point with coordinates [0,<j>). This state has energy per 
spin E = — (cos 2 + g sin 9 cos <f> + h cos 0) . States with 
spins polarised in the three orthogonal directions, |-X"+), 
\Y+) and \Z+), behave very differently with respect to 
thermalization regarding the distance between the re- 
duced evolved density matrix and the thermal one. We 
may additionally check that their dynamics are essen- 
tially different by looking at the individual expectation 
values (see also Fig. [21]). For the initial state \Y+), show- 
ing strong thermalization, all of the individual expecta- 
tion values converge fast to the thermal ones. Instead, 
for the initial state \Z+), we observe irregular oscilla- 
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FIG. 15: Distance between the averaged 3-body reduced den- 
sity matrix and the thermal state as a function of time for 
various product initial states between \X+) (/? = —0.7180) 
and \Y+) (/3 = 0).The single spin states are shown on the 
Bloch sphere on the left. 
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FIG. 17: Distance between the averaged N = 3 reduced den- 
sity matrix and the thermal one for states with similar en- 
ergy. For — 0, we compare |Y+) (red upwards-pointing 
triangles) and (9 — 0.4398,0 = n) (blue downwards-pointing 
triangles). For f3 — —0.7180, we compare \X+) (red x) to 
(9 = 1.90,0 = 0.031) (blue crosses). We show also the be- 
haviour of the state with maximal f3 — 1.5859 (green stars), 
whose expectation values oscillate like for \Z+), but in aver- 
age converges very fast to thermal. All data correspond to 
D = 120 simulation. 



FIG. 16: Distance between the averaged 3-body reduced den- 
sity matrix and the thermal state as a function of time for 
various product initial states, indicated on the Bloch sphere 
on the left, and corresponding to inverse temperatures be- 
tween j3 = -0.7180 and /3 = 0.7275. 



tions, showing no sign of damping, to the longest times 
we are able to simulate. For the initial \X+) state, we 
check that only few expectation values deviate from the 
thermal average, and are those preventing thermalization 
of the whole reduced density matrix. In particular, for 
N = 1, only (a x ) is responsible for the lack of thermal- 
ization. The behavior for larger reduced density matrices 
N = 2, 3 is qualitatively similar, although the time it 
takes for \Y+) to thermalize becomes longer. 

From the iTEBD simulations we may analyze how the 
entropy of the half-chain increases in time for the non- 
thermalizing states (see Fig. [14"]) . We notice that for the 
weak thcrmalizing initial state \Z+) the entropy of the 
half-chain grows linearly with time, so that the oscilla- 
tory behavior is not due to the absence of propagating 
excitations. For the initial state \X+) the entropy also 
grows linearly with time, but it does so at a faster pace. 
As shown in the figure, the iTEBD simulation can re- 
produce this growth until the time when the truncation 
error becomes dominant. A closer look at the data of 
this simulation shows that the distribution of Schmidt 
coefficients is very different for both cases, even at times 



when they attain a comparable entropy. If we analyze 
the Schmidt decompositions of both states at the time at 
which their entropy is S « 0.8 (t = 1.1875 for \X+) and 
t = 8 for \Z+)), we observe that the coefficient distri- 
bution for \Z+) has a much longer tail. This difference 
in the distribution can also be quantified by the 2-Rcnyi 
entropy 52 = — log(trp 2 ), which attains a lower value 
for the weak thermalizing state, S2(Z+) = 0.348, while 
S 2 (X+) = 0.611. 

By rotating the inital state on the Bloch sphere, we 
observe a transition from the strong thermalizing \Y+) 
to the weak one \Z+), and also to the apparently non- 
thermalizing \X+) (Fig. fT5]). the latter occurring around 
<t> e [f , f ] (/3 e [-0.5382, -0.3915] ). We have also anal- 
ysed the transition from weak thermalization \Z+) to 
non-thermalization \X+) (Fig. \W\i . In this case we find 
some intermediate states for which strong thermalization 
occurs. By looking at the energy, we may infer the cor- 
responding p of every initial product state. We discover 
that all the strong thcrmalizing states we have found have 
energies, and thus /3, close to zero. 

Instead, for larger (3 > 0, we observe oscillations and 
weak thermalization as in the \Z+) initial state. An in- 
teresting case is the state of maximum (3 (the product 
state with minimal energy), which we can identify, by 
studying the energy landscape over the Bloch sphere, at 
6 w 0.43, <fi = it. The dynamics of this initial state shows 
also weak thermalization, with strong oscillations of all 
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FIG. 18: Comparison between the integrable (g — 1.05, h = 
0) and non-integrable (g = —1.05, h, = 0.5) models. The plots 
show the distance between the evolved N — 3 reduced density 
matrix and the corresponding thermal one for the initial states 
\X+) (upper left), |F+) (upper right) and \Z+) (lower plot). 
The integrable case is represented always by black dots. In 
the three cases, for the integrable limit we see the relaxation 
of the reduced density matrix to a state that is different from 
thermal. We notice also that the initial state \X+) shows in 
the non-integrable case a similar behaviour to the h = limit. 



the expectation values, but fast thermalization in aver- 
age (Fig. ll7j) . We may perform a similar analysis on some 
extra states, with initial energy densities close to \X+) 
and |^+), respectively, to check whether they relax to 
similar thermal states. We find that they seem to show 
the same regime of thermalization (weak or strong) as the 
original states (Fig. IT7|) . suggesting this has to do with 
the initial energy. The relaxation curves however differ, 
indicating that, even starting with the same /3, the state 
of the system does not relax to the same state, at least 
during the long range of times we simulate. 



FIG. 19: Distance between instantaneous reduced density ma- 
trix and thermal state for N = 3, and different strengths of 
the Hamiltonian parameter g. The initial states vary from 
\Y+) to | Z+), as depicted on the Bloch sphere. Keeping 
h = 0.5 constant, we compare the case g = —1.05 (upper 
right pane) to g — —0.5 (lower left) and g — —1.5 (lower 
right). 
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Varying the Hamiltonian parameters 

Finally, we have also studied how the Hamiltonian pa- 
rameters affect the appearance of the non-thermalizing 
behaviour, to ensure that this behaviour is not singular 
to our particular choice. 

As a reference, we may compare the behavior of the 
same initial states under the chosen Hamiltonian and the 
integrable one, corresponding to h = 0. In Fig. [TS] we 
compare the dynamics under both models, h = and 
h = 0.5, for and the three most representative cases. We 
observe that in the integrable case, the N = 3 reduced 
density matrix appears to relax fast to a state which is 
not the thermal one, since the distance converges to a 
value different from zero (This could be compatible with 



FIG. 20: Distance between instantaneous reduced density ma- 
trix and thermal state for iV = 3, and different strengths of 
the Hamiltonian parameter g. The initial states vary from 
\Y+) to \X+), as depicted on the Bloch sphere. For constant 
parallel field h — 0.5, we compare the case g = —1.05 (up- 
per right pane) to g — —0.5 (lower left) and g — —1.5 (lower 
right). 



the generalized thermal ensemble [3j). 

We test also other values of the parameters g and h in 
the non-integrable regime, to study how the strong and 
weak thermalizing regimes appear. Keeping h = 0.5 con- 
stant, we observe (Fig. [T9|) that both regimes are present 
for a large range of values of g, but as we decrease g, weak 
thermalization becomes dominant, while for higher val- 
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FIG. 21: Distance between the time dependent expectation 
values of the one-body observables (a x ) (blue), (a y ) (green), 
{a z } (black) and their thermal values for the initial state \Z+) 
(left column) and \X+) (right column). Each row corresponds 
to a value of the transverse magnetic field g — —0.5 (upper- 
most), g = —1.05 (center) and g — —1.5 (bottommost). 



ues of g, the behaviour approaches strong thermalization 
in most cases. 



If we do the same now for the transition to non- 
thermalizing as observed between \Y+) and |^"+), we 
also observe that both types of behaviour survive over a 
wide range of values of g (see Fig. [20]) . 



To get a more detailed idea of the differences among the 
various types of thermalization behavior, we study the in- 
dividual expectation values for different observables. For 
clarity, we show the plots only for the N = 1 operators 
in Fig. [21] for constant parallel field h — 0.5 and varying 
g = —0.5, —1.05, —1.5, in two of the extreme cases, \Z+) 
and |^+). We observe that the oscillating behaviour of 
the initial state \Z+) appears clearly correlated with the 
value of g, and for big values the oscillations of the expec- 
tation values are clearly damped. The behaviour of the 
initial state \X+) is quite different, and in all cases we 
observe thermalization of some observables while others 
deviate from the thermal expectation value. 



